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ABSTRACT 


This  thesis  is  a  partial  analysis  of  the  Naval  Space  Command  statistical 
orbit  determination  algorithms.  Through  a  process  called  Differential  Correction, 
data  from  space  surveillance  radar  observation  stations  is  synthesized  with 
previously  accumulated  element  sets  to  maintain  accurate  orbital  object  position 
information.  Differential  Correction  is  a  nonlinear  least  squares  process  employing 
statistical  techniques  to  minimize  the  residual  measurement  error  thereby 
increasing  relative  position  information  accuracy.  This  study  focuses  specifically 
on  the  algorithmic  methods  of  solution  to  the  systems  of  normal  equations 
generated  by  the  Differential  Correction  process.  A  comparison  and  analysis  of  the 
present  Naval  Space  Command  method  and  the  singular  value  decomposition 
method  is  presented.  Algorithmic  constructions  are  presented  for  both  methods 
and  problematic  areas  are  highlighted.  The  principal  focus  herein  is  to  demonstrate 
the  benefit  of  singular  value  decomposition  when  attempting  to  solve  systems  of 
equations  whose  coefficient  matrices  are  dense  and  nearly  singular.  These  results 
generalize  to  commonly  employed  normal  equation  solution  algorithms  and  are 
intended  for  further  study  and  possible  incorporation  by  Naval  Space  Command  as 
part  of  future  modernization  plans 
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L  INTRODUCTION 


The  Naval  Space  Command  maintains  a  database  of  element  sets  for  more  than 
9,000  objects.  Through  a  process  called  Differential  Correction ,  data  from  space 
surveillance  observation  stations  is  synthesized  with  previously  accumulated  element  sets. 
Differential  Correction  is  a  nonlinear  least  squares  process  employing  statistical  techniques 
which  provide  for  accurate  orbit  state  estimation.  Radar  measurements  of  an  object’s 
motion  are  collected  by  dispersed  observation  stations  and  passed  to  NAVSPACECOM 
for  central  processing.  NAVSPACECOM  receives  daily  approximately  270,000  new 
observation  sets  and  uses  them  to  update  as  many  as  18,000  element  sets. 
NAVSPACECOM  performance  reports  indicate  that  approximately  98.5%  of  the  element 
sets  get  updated  without  manual  intervention  by  their  computer  software  called 
AUTODC.  The  remaining  1.5%  must  be  manually  analyzed.  Differential  correction  is  a 
highly  complex  step-wise  process  that  entails  more  than  16  separate  mathematical 
computations.  Even  with  automated  support,  this  constitutes  a  significant  work  load. 

Numerical  linear  algebra  is  at  the  core  of  the  differential  correction  process.  In 
general,  data  is  accumulated,  normal  equations  are  formed,  and  numerical  linear  algebra 
routines  are  called  to  solve  the  matrix  system  of  normal  equations.  The  method  of  normal 
equation  solution  and  its  associated  algorithms  are  the  focus  of  this  analysis.  The  purpose 
of  this  thesis  is  to  demonstrate  the  benefits  of  using  singular  value  decomposition  when 
obtaining  least  squares  solution  to  systems  of  normal  equations. 
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II.  BACKGROUND 


A.  STATISTICAL  ORBIT  DETERMINATION 

The  NAVSPACECOM  differential  correction  method  is  a  sequential  batch 
nonlinear  least  squares  statistical  process.  This  is  an  iterative  method  which  refines  the 
stored  orbital  element  sets  by  applying  state  adjustments  obtained  from  differential 
correction  of  the  current  orbital  observations.  These  state  updates  are  tested  for  fit  in  the 
least  squares  sense  and  if  acceptable  applied  to  the  previous  nominal  state  then  stored  as 
the  new  nominal  orbital  element  set.  The  fundamental  mathematical  steps  are  now 
outlined;  for  further  details,,  see  Vallado  [1]  and  Danielson  and  Canright  [2]. 

B.  DERIVATION  OF  THE  NORMAL  EQUATIONS 

These  descriptions  have  been  simplified  and  are  only  intended  to  provide  sufficient 
background  for  the  purpose  of  detailing  the  algorithmic  composition  and  term-wise 
structure  of  the  normal  equation  system.  For  a  very  detailed  description  of  the  entire  set 
of  NAVSPACECOM  procedures  refer  to  Danielson  and  Canright  [2]. 

C.  GENERAL  LINEAR  LEAST  SQUARES 

We  begin  the  process  by  noting  that  our  goal  is  to  solve  the  generally  inconsistent 
system  of  equations  whose  matrix  representation  is 

Ax  =b  (1) 
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by  minimizing  the  residual,  r,  where 


r  =  b  -  Ax  (2) 

The  method  of  least  squares  can  be  applied  to  solve  such  linear  systems.  Typically 
these  systems  are  highly  overdetermined.  They  have  many  more  rows  (m),  or  equations 
than  columns  (n)  or  unknowns  and  as  such  are  inconsistent.  From  elementary  linear 
algebra  [3],  we  know  that  the  range  of  the  matrix  R( A)  is  the  orthogonal  complement  of 
the  null  space  of  its  transpose  N(AT).  As  such,  the  multiplication  AT  always  produces  a 
consistent  set  of  n  equations  and  n  unknowns.  The  only  issue  with  this  method  however 
arises  when  some  set  of  the  variables  from  the  original  linear  system  are  correlated.  If  this 
happens,  the  newly  formed  consistent  system  will  fail  to  have  a  unique  solution. 

The  general  method  proceeds  as  follows.  Left  multiplication  of  equation  (1)  yields, 

AxAx  =  ATb  (3) 


Equation  (2)  implies 


Arr  =  Ax(b  -  Ax)  =  0 


Equations  from  the  original  system  may  be  weighted  by  applying  a  diagonal  weighting 
matrix.  This  expands  as 


AtWAx  =  ATWb  (5) 

Residuals  from  the  new,  weighted  system  satisfy 
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ATWr  =  ATW(b  -  Ax)  =  0 


(6) 


Provided  ArWA  is  nonsingular,  we  obtain  the  solution  as, 

x  =  (ArWA)"'  ArWb  (7) 

This  result  commonly  referred  to  as  the  normal  equations  forms  the  basis  of  the  least 
squares  solution  process. 

Differential  Correction  employs  sequential  batch  techniques  when  acquiring  the 
current  state  solution  from  different  observational  data  input  sets.  This  method 
synthesizes  the  separate  batch  solutions  as  follows. 

Given  the  systems  for  two  batches  as 

A,x  =  b,  and  A2x  =  b2 

In  order  for  the  single  solution  to  be  meaningful,  these  batch  systems  must  be 
simultaneously  solved.  This  is  accomplished  by  forming  the  weighted  least  squares  system 
as  follows. 

[(ArWA)1+(ArWA)2]x  =  (ArWb),+(ArWb)2 
The  process  can  be  simplified  slightly  by  reusing  the  solution  from  x, 

[(ArWA)j  +  (ArWA)2]x  =  (ArWA),x1  +(ArWb)2 
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When  acquiring  the  weighted  solution  simultaneously,  different  weights  may  be  applied  to 
the  different  sets  of  equations  either  by  reducing  W,  or  by  scaling  (ATWA), .  Notice  here 
that  the  batching  process  does  not  significantly  alter  the  basic  least  squares  problem  nor 
does  it  alter  the  complication  arising  when  some  portion  of  the  left  hand  side  constitutes  a 
singular  matrix. 

D  LEAST  SQUARES  APPLIED  TO  OBITAL  MECHANICS 

The  reader  is  directed  to  Danielson  and  Canright  [2]  for  extensive  mathematical 
description  and  existing  code  documentation  from  which  the  following  outline  extract  was 
taken. 


(i)  Assume  an  initial  nominal  state 


nominal 


(ii)  Compute  the  values  of  the  observed  parameters  Yc  at  N  times  corresponding  to  the 
observations  Y0  (each  observation  set  contains  at  a  maximum  6  numbers) 


,  where  /  =  1  ,...N 


(iii)  Compute  the  residuals  or  “O-Cs”  ( Ya  -  7C);  (observed  minus  calculated 
parameters).  Arrange  these  as  the  6Nxl  column  matrix 
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(iv)  Calculate  the  partial  derivatives 


SY,  37,  g(r,v) 
ax  3(r,v)  ax 

dr)  dt\ 

Hx"'~dX% 

drj  drj 

\  dYc  dYc  dYc  dYc  dYc  6YC  1  3X,  dXs 

cl  ci  c\  ci  ci  ci 

drj  drj  drK  dvj  dvj  dvK  &k 

. dX^'dX, 

. dvj  dvj 

dY^  dY^  dY^  dY^  d^_  dY^  dX,"'dXs 
drj  drj  drK  dvj  dvj  dv^J  fo,  dvj 

dx"'dXs 

wK  dvK 
_8X,  dXs_ 

at  each  observation  time  and  arrange  the  coordinates  of  the  position  and  velocity  vectors, 
(n,rj,rK)  and  (vj,vJyvK),  into  the  following  6Nx8  matrix. 


(v)  Form  the  normal  equations: 


ATWAx  =  ATWb 


Where  W  denotes  a  6Nx6N  weighting  matrix. 


Finally 


are  the  differential  corrections.  Note  that  ATWA  is  an  8x8  matrix,  and  that  ATWb  is  an 
8x1  matrix. 

(vi)  Solve  the  normal  equations: 

x  =  (Ar  WA)"1  ArWb 

(vii)  Update  the  elements: 

^new  ^old  ^ 

(For  the  first  iteration,  X0u  is  Xn0minai.) 

(viii)  Lastly,  apply  the  following  RMS  test  to  determine  if  iterations  should  continue. 

RMS,  -  RMSm  <  g 
RMS, 


where. 


RMSi 


Generally,  in  nonlinear  least  squares  our  goal  is  to  minimize  the  sum  of  the 
residuals  squared.  It  is  therefore  appropriate  to  use  some  measure  of  this  as  stopping 
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criteria.  In  the  differential  correction  process,  simply  cease  iteration  when  the  RMS 
“stops”  changing.  Since  all  of  the  necessaiy  input  to  the  RMS  equations  has  been 
accumulated  this  proves  an  efficient  method  as  well. 
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III.  NORMAL  EQUATION  SOLUTION  METHODS 

A.  GENERAL 

At  present  the  Naval  Space  Command,  Differential  Correction  process  relies  on 
Gauss-Jordan  elimination  with  full  pivoting  for  solution  to  the  normal  equations.  The 
subroutines  responsible  for  processing  the  complete  normal  equation  solutions  are  named 
AMA06,  AMA04  and  AMA03.  After  the  appropriate  Differential  Correction  algorithms 
have  formed  the  standard  system  of  normal  equations  as  described  in  the  previous  section, 
AMA06  is  called  to  read  in  the  initial  array  entries,  form  the  matrix  E  =  ArWA ,  duplicate 
it  and  begin  preconditioning  processes  designed  to  identify  singular  matrices.  Following 
sufficient  preconditioning,  AMA06  then  calls  AMA04  which  in  turn  computes  E”1 . 

Finally,  AMA03  is  called  to  apply  Gauss-Jordan  elimination  and  solve, 
x  =  (ArWA)-1  ArWb .  The  8x8  matrix  ArWA  is  input  as  E  and  the  8xlmatrix  ArWb 
is  input  as  G. 

The  following  section  describes  specific  processes  and  algorithms  used  to  solve  the 
normal  equations  with  Gaussian  elimination  and  with  Singular  Value  Decomposition, 
SVD.  Each  process  is  initiated  with  the  previous  normal  equation  derivation.  The  present 
NAVSPACECOM  differential  correction  process  uses  the  Gauss  Jordan  elimination 
method  with  full  pivoting. 
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Focusing  back  on  step  (vi)  of  the  least  squares  orbital  mechanics  process,  solution 
of  the  normal  equations,  we  now  compare  the  present  NAVSPACECOM  method  with 
Singular  Value  Decomposition  method.  The  refined  problem  is  now,  how  to  acquire  the 
best  approximation  at  each  step  of  the  nonlinear  approximation  process  so  as  to  ensure 
maximum  efficiency  of  the  iterative  routine  while  minimizing  the  associated  rounding 
error  in  the  next  iterate? 

B.  NAVSPACECOM  GAUSS-JORDAN  ELIMINATION 

The  Naval  Space  Command  uses  Gauss-Jordan  elimination  with  full  pivoting  for 
solution  to  the  system  of  normal  equations.  Their  algorithm  is  very  similar  to  the  one  in 
the  book,  Numerical  Recipes,  by  Press,  et  al.  [7], 

Typical  Gauss-Jordan  with  Full  Pivoting 

Eveiy  Gauss-Jordan  step  is  a  left  multiplication  by  an  elementary  matrix,  [3],  An 
overview  of  the  general  method  of  Gauss-Jordan  elimination  is  as  follows: 

(i)  Augment  the  right  hand  side  with  the  nxn  identity  matrix. 

(ii)  Perform  elementary  operations  on  the  augmented,  nx2n  system  until  the 
nxn  identity  matrix  is  reformed  in  the  first  n  columns  of  the  matrix 

(iii)  Recall  the  last  n  columns  of  the  augmented  matrix  as  the  inverse. 
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This  method  is  usually  done  in  place;  that  is,  the  actual  augmentation  is  never 
performed.  The  algorithm  simply  replaces  the  original  input  matrix,  in  place,  with  its 
inverse.  For  this  reason,  if  there  is  ever  a  need  to  recall  the  original  data,  a  copy  of  the 
original  matrix  must  be  made  on  input. 

Full  pivoting  is  the  process  by  which  rows  and  columns  of  the  original  matrix  are 
reordered  so  that  the  element  largest  in  magnitude  is  moved  to  the  upper  left  comer  of  the 
matrix.  Then,  each  successive  pivot  row  in  maneuvered  similarly.  Why  do  this?  In  the 
computational  process,  the  subsequent  rows  below  the  pivot  in  question  are  being  reduced 
to  zero  by  dividing  through  by  a  scaled  multiple  of  that  pivot.  If  the  process  is  not  begun 
in  this  manner,  that  is  to  say  if  the  element  below  the  pivot  is  greater  in  magnitude  then  the 
pivot,  the  pivot  must  be  scaled  by  increasing  its  magnitude.  The  scaling  factor  is  the 
element  being  zeroed  out  divided  by  the  pivot,  and  if  the  pivot  is  small  enough  this 
approaches  division  by  zero.  In  finite  precision  arithmetic,  this  can  lead  to  a  divide  by  zero 
condition.  If  the  difference  in  magnitude  of  the  pivot  is  sufficiently  less  than  the  element 
being  zeroed  out ,  the  computer  can  evaluate  the  expression  and  return  “Not-A-Number”. 
Even  if  the  division  by  zero  extreme  condition  does  not  occur,  scaling  the  matrix  through 
by  these  sufficiently  small  quantities  inflates  the  rounding  error.  How  does  pivoting 
protect  against  this?  By  always  scaling  so  that  each  successive  pivot  is  the  next  largest  in 
magnitude,  the  previous  process  of  dividing  the  subdiagonal  elements,  results  in  scaling 
the  pivot  row  by  decreasing  its  magnitude.  This  has  the  effect  of  moving  the  division 
operations  away  form  the  1/zero  condition  and  deflates  the  magnitude  of  rounding  error  to 
the  minimum  possible  given  the  actual  element-wise  composition  of  the  matrix. 
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Partial  pivoting  is  a  similar  process  that  allows  only  row  exchanges.  While  this 
may  appear  to  be  a  sound  concept  with  less  computational  overhead,  in  practice,  matrix 
elements,  not  directly  below  (in  some  other  column)  the  pivot  being  operated  on,  may  still 
be  sufficiently  greater  in  value  than  subsequently  available  pivot  choices.  If  this  occurs, 
the  solution  is  again  driven  toward  a  divide  by  zero  condition.  By  manipulating  both  rows 
and  columns,  as  in  foil  pivoting,  we  reduce  the  likelihood  of  the  divide  by  zero  condition 
to  the  greatest  extent  possible. 

When  employing  the  full  pivoting  routines,  the  general  Gauss-Jordan  process  is 
augmented  by  left  and  right  permutation  matrices  which  track  the  pivot  maneuvers.  The 
corresponding  information  is  stored  as  an  array  and  applied  during  the  back  solve  process 
thereby  ensuring  correct  association  of  the  coefficients  and  variables  from  the  original 
system  of  equations. 

It  must  be  noted  that  the  previous  pivoting  processes  only  reduce  the  likelihood  of 
algorithm  failure.  They  do  not  prevent  it.  In  the  event  of  near  singular,  noninvertible 
conditions,  the  Gauss-Jordan  full  pivoting  routine  has  no  mechanism  to  correct  or 
compensate.  This  is  an  inherent  critical  deficiency.  Failure  due  to  nearly  singular  input 
matrices  routinely  occurs  when  processing  least  squares  systems  that  are  derived  from 
closely  sampled  data. 
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Unique  features  of  NAVSPACECOM  Gauss-Jordan 


The  following  system  processes  are  extracted  as  reported  in  the  code 
documentation  research  of  Danielson  and  Canright,  [2],  The  matrix  representation  of  the 
system  of  normal  equations  is  loaded  as  the  input  matrix  E.  The  input  matrix  is  symmetric 
positive  semi-definite  by  its  inherent  construction.  As  such,  only  values  on  and  to  the  right 
of  the  diagonal  are  read  as  input.  The  below  diagonal  elements  are  copied  from  their 
respective  symmetric  counterparts.  At  this  point  a  preconditioning  routine,  AMA06,  is 
called.  It  samples  the  matrix  to  determine  if  any  off-diagonal  element  of  E  is  too  large 
relative  to  its  corresponding  diagonal  elements.  AMA06  tests  each  E?  >  (SJNG)EnEjl , 

where  SING  is  a  parameter  value.  If  any  E?  violates  this  inequality  then  that  row  (for 

even-number  calls)  or  column  (for  odd-number  calls)  is  “ inactivated The  threshold 
SING  is  initially  set  to  near  one.  Then,  if  the  active  matrix  is  singular,  the  threshold  SING 
is  lowered  by  10%  (attempting  to  inactivate  more  rows/columns)  and  the  solution  is  tried 
again.  At  this  point,  a  saved  copy  of  the  original  matrix  must  be  reloaded.  The  immediate 
effect  of  this  process  is  to  eliminate  divide  by  zero  conditions  when  pivot  elements  are 
very  small.  Eventually,  either  a  solution  is  found,  or  the  whole  matrix  is  made  inactive 
(flagged  by  a  return  value  SING=0),  or  the  threshold  gets  too  small  (SING  <0.01) 
indicating  that  the  entire  system  is  singular.  In  the  event  that  no  matrix  remains,  an  error 
condition  is  returned  and  all  further  processing  of  this  set  of  observational  data  ceases.  A 
tally  of  these  failures  is  presented  as  an  output  report  in  the  batch  run  statistics. 
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There  is  a  fundamental  problem  with  this  preconditioning  methodology.  Under  no 
circumstances  are  the  root  causes  of  the  perceived  singularities  dealt  with.  Further,  in 
terms  of  computational  effort,  the  computer  is  permitted  to  cycle  excessively  without 
sufficient  intermediate  checks  to  determine  if  the  initial  perceived  singular  conditions  are 
removable. 

In  summary,  Gaussian  elimination  is  generally  the  most  computationally  efficient 
numerical  method  available;  however,  it  suffers' from  numerical  instability  under  certain 
conditions,  which  are  detailed  in  Chapter  IV,  Comparison  of  Methods.  Unfortunately, 
these  very  same  conditions  generally  arise  in  systems  of  normal  equations.  As  noted 
earlier,  whenever  correlation  of  variables  from  the  original  system  equations  occurs,  the 
resulting  matrix  Ar  A  represents  a  consistent  set  of  equations  with  redundant  solutions. 
These  redundant  solutions  drive  the  matrix  toward  a  singular  condition. 

C.  SINGULAR  VALUE  DECOMPOSITION 

The  Singular  Value  Decomposition,  SVD  of  an  arbitrary  mxn  matrix  is  the 
factorization  of  A  into  UIVT,  where  U  and  V  are  orthogonal  matrices  and  Z  is  the 
rectangular  mxn  matrix  whose  first  r  rows  form  a  square,  diagonal  submatrix  with 
elements  ax---ar  ,  i.e.,  the  singular  elements  of  A  with  the  remainder  of  Z  being  zeros. 

The  process  proceeds  generally  as  follows  (See  [6]): 

(i)  Reduce  the  general  input  matrix  A  to  bidiagonal  form  B  with  orthogonal 

matrices  U  and  V  where  A  =  U  BVr  .  B  is  nonzero  only  on  its  main 
11  11 
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diagonal  and  first  super  diagonal.  This  is  accomplished  with  Householder 
transformations  in  the  actual  algorithms. 

T 

(ii)  Find  the  SVD  of  B,  B  =  U  £V  ,  where  £  is  the  diagonal  matrix  of  singular 

2  2  W 

values  and  and  V  are  orthogonal  matrices  whose  columns  are  the 

respective  left  and  right  singular  vectors.  The  Gram-Schmidt  process 
produces  this  result. 

(iii)  Combine  these  decompositions  to  form  A  =  (U^Up^VV^/ .  The 
columns  of  U  =  (UU2)and  V  =  (V  V  )are  the  respective  left  and  right 
singular  vectors  of  A. 

Step  (i)  reduces  the  A  matrix  to  bidiagonal  form  by  applying  Householder 
transformations  on  both  left  and  right  sides.  The  symmetry  of  the  original  matrix  is 
preserved  throughout.  Step  (ii)  employs  the  Gram-Schmidt  process  to  orthogonalize  the 
columns  of  and  V  .  Step  (iii)  reforms  the  matrix  factors.  Specific  algorithms  for 

these  subordinate  functions  are  detailed  in  Golub  [5]  and  Demmel  [6],  The  SVD 
algorithm  simply  incorporates  them.  Another  extremely  useful  secondary  benefit  of  the 
SVD  is  provided  by  its  fundamental  structure.  That  is  the  factor  matrices  U  and  V  have  a 
very  special  structure.  They  are  orthogonal.  Where,  A  =  U£VT  ,  the  columns  of  U  are 
the  eigenvectors  of  AAT  and  the  columns  of  V  are  the  eigenvectors  of  ATA  ,  Strang, 

[3],  Further,  orthogonal  matrices  have  the  nicety  that  U  T  U  =  UUT  =1.  This  infers 
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directly  that  when  solving  the  system  of  normal  equations. 

Ax  =  b  as  UEVtx  =  b  the  respective  transposes,  UT  and  V,  are  left 
multiplied  yielding  x  =  VI+UTb .  Interestingly  enough,  no  actual  inverse  matrices  have  to 
be  calculated.  Further,  the  transpose  of  U  and  V  can  be  done  in  place.  The  real  benefit  is 
the  complete  absence  of  rounding  error  when  taking  the  transpose  of  a  matrix. 

The  problem  with  the  Gauss-Jordan  method  is  the  requirement  for  the  columns  of 
A  to  be  independent.  As  described  in  detail  in  the  previous  section,  when  ill-conditioned 
matrices  are  input,  a  great  deal  of  computational  effort  is  required  in  preconditioning  the 
matrix.  Without  this  preconditioning,  even  nonsingular  but  very  ill-conditioned  matrices 
may  cause  algorithms  to  break  down.  The  computer  simply  perceives  the  matrix  to  be 
singular  to  its  level  of  machine  precision. 

The  key  to  SVD’s  computational  stability  lies  in  its  orthogonality.  Matrix  rank 
problems  arise  frequently  in  computer  arithmetic.  Determining  the  rank  of  an  ill- 
conditioned  matrix  can  be  challenging  in  the  presence  of  roundoff  error  and  noisy  data. 
The  SVD  allows  for  practical  dealing  with  numerical  rank  deficiency 

The  following  theorem  is  taken  from  Golub,  [5],  See  [5]  for  a  detailed  proof.  The 
theorem  provides  the  detail  necessary  for  determining  how  “close”  the  given  A  matrix  is  to 
one  of  lower  rank.  The  2-norm  here  is  the  matrix  derivation  of  the  vector  Euclidean  norm 
[3]  defined  as  follows. 


least  upper  bound 


(A)=max 


^"Axl1^ 


V  Fli  J 


■ where 


A  <  least  upper  bound (A)||x[[ 
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Theorem 


Let  the  SVD  of  A  eRm“  be  given.  Ifk  <  r  -  rank(A),  and 

k  ' 

Ak  = 

<=i 

then, 

min  ||a-B||  =  A- A  =crk  +  ] 

11  "2  k  ||  2 

rank(B)  =  k 

This  striking  result  offers  a  method  of  computation  using  stored  values  which 
specifically  reveal  the  magnitude  of  the  degenerate  numerically  singular  condition  of  the  A 
matrix.  The  details  of  this  theorem  indicate  that  the  smallest  singular  value  of  A  is  the  2- 
norm  distance  of  A  to  the  set  of  all  rank  deficient  matrices.  Iteration  is  no  longer  required 
to  determine  the  degree  of  singularity. 

Now,  let’s  look  at  an  example  of  how  to  make  the  most  out  of  ill-conditioned  least 
squares  systems  with  the  SVD.  See  Strang  [3]  for  more  on  this.  In  general,  the  least 
squares  problem  has  one  very  stringent  requirement,  the  columns  of  the  A  matrix  must  be 
independent  or  the  rank  of  A  must  be  equal  to  n,  the  number  of  columns.  This  is  often 
referred  to  as  full  rank.  If  not,  A  is  not  invertible  then  Ax  =  b  can  not  determine  x.  As 
described  earlier,  any  vector  from  the  null  space  of  A  can  be  added  to  x.  Now  let’s 
examine  what  happens.  There  are  two  possible  situations,  either  the  rows  of  A  may  be 
dependent  or  the  columns  of  A  may  be  dependent.  The  first  situation  implies  the  system 
of  equations  may  have  no  solution  and  the  second  situation  implies  that  any  solution  is  not 
unique.  The  dependent  column  case  makes  this  a  particularly  difficult  yet  interesting 
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problem.  As  discussed  earlier,  when  we  have  dependent  rows  the  solution  we  seek  may 
be  outside  the  column  space  of  A.  Our  course  of  action  now  becomes  simply  project  b 
onto  the  column  space  of  A.  Now  the  greater  challenge.  After  making  that  projection  we 
find  A  has  dependent  columns  and  the  solution  is  not  unique.  At  this  point,  we  must  now 
employ  the  criteria  for  selecting  the  optimal  solution  and  choose  the  one  with  minimum 
length. 

Consider  the  following  example: 

cr,  0  0  0 

0  cr2  0  0 

0  0  0  0 

0  0  0  0_ 

where  A  is  diagonal  with  dependent  rows  and  columns. 

Here  we  see  the  columns  all  end  in  0.  As  per  case  (i),  the  closest  vector  to 
b  =  (&i,Z>2,Z>3,Z>4  )  is  p  =  0,0) the  projection  onto  the  column  space  of  A.  The 

magnitude  of  error  here  is  b  =  (0,0,Z>3,Z>4)the  perpendicular  to  the  columns.  The  best 
solution  now  is  attained  when  we  solve  the  first  two  equations.  Since  the  last  two 
equations  indicate  0  =  b3  and  0  =  Z>4,  the  error  in  those  equations  cannot  be  reduced  but 
the  error  in  the  first  two  will  be  zero. 


Ax  =  p  is 


cr1  0  0  0 

0  <r2  0  0 

0  0  0  0 

0  0  0  0 


JCi 

bx 

X2 

b2 

Xi 

0 

X4_ 

0_ 

(8) 
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Now  the  challenge,  the  dependent  columns  imply  x  is  not  unique!  The  first  two 
components  are  —  and  — ,  but  x3  andx4  are  completely  arbitrary.  Now  apply  the 

minimum  length  criteria  and  see  that  these  arbitrary  components  must  be  identically  zero 
to  attain  the  best  approximation. 


That  is. 


(9) 


The  minimum  length  solution  to  Ax  =  p  is ,  x+  ,  Strang  [3],  Again,  the  useful  result  is 
specifically  the  equation  that  reveals  x+  .  This  process  displays  the  matrix,  which  yields 
the  desired  result, 


//A 


0 

0 

0 


0 

^2 

0 

0 


0  o' 

0  0 

then  A 

0  0 
0  0_ 


—  000 
0i 

0  —  0  0  and  x+ 

0  0  0  0 

0  0  0  0 


(10) 
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A+  is  referred  to  as  the  pseudoinverse  and  is  the  matrix  which  provides  for  solution  to  the 
nearly  singular  system  of  Ax  =  b  ,  Strang,  [3], 

This  entire  process  has  one  sticking  point.  £+  is  the  pseudoinverse  described 
earlier.  Now,  the  magnitude  of  rounding  error  in  applying  the  inverse  process  is  limited  to 
the  sum  of  the  errors  when  inverting  the  individual  diagonal  elements,  cr. .  Each  <7i  is  the 

square  root  of  each  nonzero  eigenvalue  from  both  AAT  and  ATA  .  Now,  the  fine 
point,  what  happens  when  cr.  is  sufficiently  small  to  induce  a  divide  by  zero  condition 
when  taking  the  pseudoinverse?  The  reciprocal  of  cr,  is  set  to  to  zero  by  the  code.  Press, 

et.  al.,  [7],  denote  this  procedure  as  editing  the  singular  values.  The  logic  is  sound. 

Recall  the  original  formulation  of  the  linear  system.  When  redundant  solutions  (singular 
conditions)  are  encountered  as  a  result  of  variable  correlation,  the  matrix  is  unable  to 
distinguish  between  the  different  basis  functions  and  the  associated  distribution  of  the 
input  data.  By  setting  the  reciprocal  of  any  sufficiently  close  to  zero  singular  value  to 
zero,  we  effectively  add  a  zero  multiple  to  the  fitting  parameters  as  opposed  to  some  large 
combination  of  the  basis  functions  that  are  degenerate  to  the  best  fit,  [7],  Further,  if  any 
nonzero  singular  value  is  very  small,  its  reciprocal  should  also  be  set  to  zero.  This  term  is 
most  likely  residual  from  rounding  error  and  detracts  from  the  optimal  solution. 

A  rule  for  determining  the  editing  tolerance  of  singular  values  is  given  by  Press, 
et.  al.,  [7]  as  follows.  Set  to  zero  any  singular  value,  cr,.  when, 
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The  ratio 


<  Ns 


where  N  is  the  column  length  dimension  and  s  is  the  machine  precision. 

Consider  the  following  example  least  squares  problem  which  illustrates  the 
mathematical  principles  of  the  algorithm.  In  this  example  A  is  assumed  symmetric  on 
input  just  as  in  the  Differential  Correction  form. 


where 


Let  A  = 


3-2  2 
-2  4  0 

2 


and  b= 


A^A  = 


0 

2J 

Li 

'  17 

-14 

10" 

-14 

20 

-4 

10 

-4 

8 

The  eigenvalues 


X  =  36=6  =><j  =6 

1  i 

A^A  -  XI  are  X  =9=3 2  =>cr  =3 

2  2 

1  =0=02^C7  =0 

3  3 


whose  corresponding  eigenvectors  are 


1 

v  =  — 

'  2  " 

-2 

1 

,  v  =  — 

T 

2 

1 

,  v  =- 

'  2  ' 

1 

1  3 

1 

2  3 

2 

3  3 

-2 

Now  A  =  UZVror  AV=UE  so 
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u  = — Av  ,  u  = — Av 
i  cr  1  2  cr  2 

1  2 


'3  -2  2 

-2  4  0 

'  2  “ 

-2 

1  _  1 

'12' 

12 

_  1 

«  7 

_  2  0  2_ 

_  1  _ 

3  “  18 

_  6  _ 

”  3 

_  1  _ 

Similarly 


"3  -2  2" 

-2  4  0 

T 

2 

1  _  1 

'3' 

3 

1 

T 

2 

1 - 

to 

o 

. 

2_ 

3  "  9 

6_ 

3 

2_ 

Now  for  the  zero  singular  value  we  must  solve  the  homogeneous  problem 
A  ru  =  0  for  . 

-2  2 


'3  -2  2 

-2  4  0 

o  o 

3  3 

0  1  - 

2  0  2 

0 

2 

0  0  0 

1 


0 

0 

0 


is  arbitrary 


but ,  u  =—u  and  u  =-u 

2  2  3  13 


so. 


u  -un 
3  3 


’-l" 

'  2  " 

1 

_  1 

1 

2 

"  3 

1  - 

_-  2_ 

Now  lets  refit  the  components. 
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r  ,[6  °1  v> 

A  =  U£Vr  =  u  u  u  0  3  0  vr 

1  2  3J  2 

0  0  OJ  vr 

|_  3  _ 

This  is  simply  the  sum  of  rank  one  matrices 
Now, 

2  -1 

3  3 

-2  -2 

3  3 

I  z2l 

3  3 

Now  for  the  least  squares  piece. 

Ax=b  and  A  =  UEVr  so  2Vrx  =  Urb  or  x  =  VS_1Urb 

IT 1  is  the  pseuodoinverse.  Its  structure  was  detailed  earlier. 

Now  applying  equation  (8)  to  the  original  right  hand  side  yields, 

>  0  I]  [2] 

9  9  1  9 

0  -  -  1  =  - 

9  9  3 

1  I  I L1-!  J_ 

9  9  6J  [l8. 

This  is  the  least  squares  solution  of  minimal  length. 

Now  the  test,  if  x  is  orthogonal  to  N(  Ar ),  i.e.,  really  the  minimal  length  solution 
and  u3  is  an  element  of  N(AT  )  then  orthogonality  demands  that  their  dot  product  equal 
zero. 


A  =  U2Vr=> 
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This  example  is  easily  extended  to  matrices  of  greater  dimension. 
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IV.  COMPARISON  OF  METHODS 


A.  GENERAL 

The  procedures  used  for  computational  comparison  were  NAVSPACECOMs 
algorithms  AMA03,  AMA04,  and  AMA06  coupled  with  a  process  driver  program  and  the 
standard  SVD  algorithms  as  taken  from  Numerical  Recipes  bv  Press,  et.  al.,  [7].  During 
research  of  the  actual  NAVSPACECOM  code,  it  was  found  that  Numerical  Recipes  was 
referenced  repeatedly  throughout.  This  observation  set  the  precedence  for  applying  the 
basic  SVD  codes  of  Numerical.  Recipes  for  computational  analysis.  This  common  ground 
should  serve  to  standardize  any  resulting  code  development. 

It  should  be  noted  that  throughout  chapter  15,  Numerical  Recipes  [7]  the  SVD  is 
recommended  for  least  squares  problems.  Citing  the  following,  “...solution  of  a  least 
squares  problem  directly  from  the  normal  equations  is  rather  susceptible  to  roundoff 
error.”  And,  ”In  some  applications,  the  normal  equations  are  perfectly  acceptable  for 
linear  least  squares  problems.  However,  in  many  cases  the  normal  equations  are  very  close 
to  singular.”  The  authors  detail  with  significant  correlation  the  precise  difficulties  that 
generally  arise  in  Differential  Correction.  They  go  on  further  detailing  the  exact 
complication  the  routine  AMA06  attempts  to  eliminate.  The  following  excerpt  from 
Numerical  Recipes  sums  up  the  exact  nature  of  this  entire  analysis: 
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A  zero  pivot  element  may  be  encountered  during  the  solution  of  the  normal 
equations,  with  Gauss-Jordan,  in  which  you  get  no  solution  at  all.  Or,  a  veiy  small  pivot 
may  occur  in  which  case  you  typically  get  fitted  parameters  ak  with  very  large 

magnitudes  that  are  delicately  balanced  to  cancel  out  almost  precisely  when  the  fitted 
function  is  evaluated. 

Why  does  this  commonly  occur?  The  reason  is  that,  data  do  not  clearly 
distinguish  between  two  or  more  basis  functions  provided.  If  two  such  functions  or  two 
different  combinations  of  functions  happen  to  fit  the  data  equally  well  -  or  badly  -  then 
the  matrix  A,  is  unable  to  distinguish  between  them  and  becomes  singular.  There  is 
irony  in  the  fact  that  least  squares  problems  are  both  overdetermined  (number  of  data 
points  greater  than  the  number  of  parameters)  and  underdetermined  (ambiguous 
combinations  of  parameters  exist).  The  ambiguities  can  be  extremely  hard  to  notice  a 
priori  in  complicated  problems. 

The  SVD  gives  exactly  what  we  need.  In  the  overdetermined  system  SVD 
produces  a  solution  that  is  the  best  approximation  in  the  least  squares  sense.  In  the  case 
of  the  underdetermined  system,  SVD  produces  a  solution  whose  values  (the  ak ’s)  are 

the  smallest  in  the  least  squares  sense  also  what  we  want.  When  some  combination  of 
the  basis  functions  is  irrelevant  to  the  fit,  that  combination  is  driven  down  to  a  small, 
innocuous  value,  rather  than  pushed  up  to  delicately  canceling  infinities. 


B.  SENSITIVITY  ANALYSIS 

The  following  sections  present  information  that  details  the  inner  workings  of  the 
individual  numerical  solvers.  The  concept  to  focus  on  deals  with  the  aspect  of  matrix 
sensitivity.  The  following  simplified  example  details  the  problem  explicitly  (see  Nyhoff 
[16]). 


Ax  =  b  is 


1 - 

VO 

t<N 

X 

8 

[2  6.0000003J 

y_ 

_8.0000003_ 

The  solution  is  obviously  x  = 


.  Now  perturb  the  right  hand  side  very  slightly. 
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Ax  =  b*  is 


'2 

6 

X 

8 

_2 

6.0000003_ 

_y_ 

_7.9999994_ 

Now  the  not  so  obvious  solution  is  x  = 

|  right  hand  side  has  altered  the  solution  on  the  order  of  107  times  the  constant  term  of  the 

! 

|  second  equation. 

In  performing  rounding  error  analysis,  we  note  that  computers  operate  in  floating 
point  arithmetic.  That  is,  they  convert  every  problem  into  a  “ nearby ”  problem  perform 
calculations  and  then  convert  back  the  answer.  The  symbol  s,  called  machine  epsilon 
refers  to  the  computer’s  ability  to  distinguish  between  two  consecutive  binary 
representations  of  actual  input  values.  When  the  relative  representation  of  two 
consecutive  numbers  exceeds  the  computer’s  ability  to  distinguish  between  them, 

(  |  a- Z>||  <  s  )  we  say  an  underflow  has  occurred  in  the  calculation  sequence. 

The  following  method,  known  as  matrix  perturbation  theory,  helps  to  analyze  the 

i 

nature  of  error  in  order  to  minimize  it  through  algorithm  refinement.  Errors  are 
accumulated  primarily  from  two  sources.  First  there  may  be  measurement  error  contained 
within  the  original  input  data,  and  second,  the  algorithm  itself  may  cause  error  through  its 
internal  approximations  while  processing  calculations.  For  excellent  derivations  of  the 
most  commonly  required  error  measurement  techniques,  see  Demmel  [6]. 


10 

-2 


,  where  the  very  small  perturbation  in  the 
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The  measure  of  the  accumulated  error  is  referred  to  as  the  condition  number  of  a 
matrix.  Simply  put,  the  condition  number  indicates  the  precise  level  of  sensitivity  a  matrix 
has  relative  to  these  types  of  very  small  perturbations. 

To  summarize,  an  example  taken  from  Strang  [3]  is  in  order.  It  directly  illustrates 
the  short  fall  in  the  present  NAVSPACECOM  code  methodology  of  scaling  the  input 
matrix  by  the  value  25,000,000. 

Begin  by  consider  the  linear  system  Ax  =  b ,  now  perturb  the  right  hand  side  by 
5  b  .  These  errors  might  have  come  from  the  observational  data  or  roundoff.  The  change 
is  small  but  the  direction  of  change  can  not  be  controlled.  The  solution  is  subsequently 
changed  from  x  to  x  +  S  x .  Now  our  system  has  become  AS  x  =  <5b .  We  now  must 

estimate  the  resulting  perturbation  Sx  =  A-1<5b .  There  will  be  a  large  change  in  the 

solution  when  A-1  is  large,  i.e.  A  is  nearly  singular.  Now  consider  our  symmetric  matrix 
with  positive  eigenvalues  where 0< /l,  <  X,  <•■■<!„.  Any  vector  <5b is  a  combination  of 
the  corresponding  unit  eigenvectors  X!,"-,xn .  Let  e  indicate  a  very  small  change. 

If  <5b  =  £x,  then  Sx  =  — . 

4 

The  error  in  ||<5b||  is  amplified  by  — ,  which  is  the  largest  eigenvalue  of  A'1 .  The 
amplification  is  greatest  when  \  is  close  to  zero.  Thus  nearly  singular  matrices  are  the 
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most  sensitive.  The  serious  drawback  with  this  measure  of  sensitvity  appears  when  scaling 
is  introduced.  Multiplying  the  matrix  by  a  large  scalar  also  scales  the  eigenvalue  by  the 
corresponding  amount.  This  makes  the  matrix  appear  much  less  singular.  This  rescaling 
can’t  however  make  an  ill-conditioned  matrix  well.  It  is  true  <5x  will  be  smaller  by  the 


scale  factor  however  so  will  the  solution  to  x  =  A  'b .  The  relative  error 


.stays  the 


same.  The  factor  ||x||  in  the  denominator  normalizes  the  problem  against  such  trivial 
rescaling.  There  is  a  corresponding  normalization  of  ||<5b|| .  The  problem  is  to  compare  the 


relative  change 


with  the  relative  error 


|<5x 


At  this  point,  I  refer  to  the  following  theorem  as  taken  from  Strang,  [3], 


Theorem 


For  a  positive  definite  matrix,  the  solution  x  =  A  :b  and  the  error  5  x  =  A  ’(5b  always 
satisfy, 


HI  *  j-  and  INH 

Therefore  the  relative  error  is  hounded  by 
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JMI<iJNI 
Ml  '  4  IMI  ' 

X  X 

The  ratio  c  -  — -  =  — 'SSL  is  called  the  condition  number  of  A. 

\  Anin 

This  analysis  can  be  applied  to  the  example  given  previously. 

C.  ERROR  ANALYSIS 

The  following  analysis  applies  to  the  least  squares  system.  The  established  error 
bounds  hold  regardless  of  solution  method.  A  fundamental  rounding  error  analysis 
follows  directly  from  the  previous  sensitivity  analysis.  See  Golub,  [5]  and  Demmel,  [6]  for 
very  complete  presentations  of  both  rounding  and  backward  error  analysis 

The  tradition  definition  of  a  vector  or  matrix  norm  holds  throughout  the  following 
analysis.  See  Strang  [3]  as  required.  Where  the  symbol  ||.  |]  appears  in  equation  form, 

consistent  use  of  the  chosen  norm  is  implied  throughout.  Condition  number  is  as  defined 
previously. 

Given  the  linear  system  Ax  =  b  where  r  =  b  -  Ax  .  Where  A  is  nxn 
nonsingular,  x  is  the  computed  solution,  and  r  is  the  residual.  Letxe  be  the  exact  solution 
and  e  the  error. 

Then  r  -  b  -  Ax  =  Axe  -  Ax  =  A(xe  -  x) 
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or 


r  -  Ae 


Given  A  nonsingular  then  e  =  A  V .  Taking  consistent  norms  of  the  equations  yields  the 
following  inequalities. 


< 


HI 


Relating  these  inequalities  bounds  the  error  as  follows, 


. -i 


r  >  e  > 


(12) 


Now  the  exact  solutions  form, 


xe  =  A  'b  and  Axe  =  b 


Where 


<  A“ 


b  and  A  kj  >  b 


We  can  now  bound  the  exact  solution. 


Now  divide  inequality  (10)  through  by  \xe  and  form 
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> 


Substituting  llxj  in  the  left  denominator  by  f- JL  and  in  the  right  denomination  by 

iAll 

( 

IIa-1  |jbj|  the  expression  for  relative  error  now  becomes, 


(13) 


Where  tt-4  is  the  relative  error. 


Now  recalling  that  the  condition  number  of  the  matrix  A  can  also  be  expressed  as 


cond( A)  =  ||a-1||||a||  ,  which  is  required  for  analysis  of  methods  that  do  not  return 
eignevalue  information,  equation  (11)  becomes. 


cond(  A) 


fr\) 

>  HI  >  i 

M 

lllbJ 

|xe||  cond  (A) 

lbJ 

This  formulation  is  convenient  in  that  it  allows  for  testing  of  the  maximum  and 
minimum  associated  errors  for  any  given  output  as  initiated  by  the  condition  of  the  matrix 
system.  For  more  specific  information  on  relative  error  and  on  absolute  error  given 


34 


machine  dependent  performance  information  see  Engeln-Mullges  and  Uhlig  [7],  This 
analysis  can  be  applied  to  the  given  test  matrices. 

D.  COMPUTATIONAL  EFFICIENCY 

Flops,  or  floating  point  operations,  are  the  measure  of  algorithm  efficiency.  While 
they  do  not  necessarily  indicate  the  actual  processing  time,  as  different  computers  do 
internal  processing  differently,  they  do  give  an  excellent  measure  of  the  magnitude  of  total 
number  of  computer  calculation  required  by  an  algorithm  as  related  to  the  functional  input. 

In  general,  flop  counts  are  obtained  by  summing  the  number  of  arithmetic 
operations  from  the  most  deeply  nested  algorithm  statements.  As  an  example  of  the 
accounting  notice,  a  dot  product  operation  of  length  n  involves  n  multiplications  and  n 
additions.  It  therefore  requires  2n  flops.  For  matrix  multiplications,  the  general  form  is 

C(/,y)  =  A(i,k)*B(k,j)  +  C(i,j) 

This  process  requires  2mnp  flops  where  C  e  R"“" ,  A  sR'”7’,  B  e  R/,x” . 

Applying  the  above  flop  counting  procedure,  the  following  flop  count  estimates 
are  given  for  each  algorithm  as  related  to  the  matrix  dimension  [5], 
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ALGORITHM 

FLOP  COUNT 

Gauss-Jordan 

2  «3 
mn  +  — 

3 

SVD 

4  m2n  +  Smn2  +  9n 3 

(U,2,V) 

SVD 

4  mn2  +  8  n2 

(2,  V  only) 

Minimum  requirement  for  least  squares 

This  does  not  take  into  account  specialized  storage  techniques  and  factorizations 
when  dealing  with  specific  subprocess  forms  that  may  be  optimized.  In  this  respect,  these 
estimates  would  be  considered  worst  case.  In  some  instances  the  total  flop  count  of  a 
particular  subform  may  be  cut  by  nearly  half  the  indicated  estimate.  While  these  methods 
yield  respectable  approximations  of  the  computational  cost,  they  do  not  provide  actual  run 
time  analysis.  Specific  algorithms  can  be  made  extremely  efficient  when  optimized  for  a 
particular  computer. 
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V.  CONCLUSIONS 


A.  SUMMARY  OF  FINDINGS 

Testing  of  the  actual  NAVSPACECOM  code  and  the  SVD  code  was 
accomplished  through  PC  based  Digital  Visual  Fortran  6.0.  All  NAVSPACECOM  source 
code  was  compiled  and  linked  using  the  visual  Fortran  development  environment.  As 
noted  earlier  SVD  source  codes  were  adopted  directly  from  Numerical  Recipes  [6]. 

Driver  programs  for  both  software  suites  were  developed  as  adoptations  of  the  LAPACK 
[9]  and  Numerical  Recipes  source  code  for  advanced  linear  algebra.  Test  matrices  were 
generated  using  MATLAB  and  placed  in  standard  text  format  for  file  upload  by  the  driver 
programs.  All  Fortran  source  code  was  version  77.  The  computer  system  used  was  a 
Micron,  300  MHz,  Pentium  PH. 

The  test  matrices  indicated  at  Appendices  A  and  B  were  used  to  determine  how 
well  the  two  routines  compared  when  attempting  to  solve  highly  singular  and  known 
singular  input  matrices.  As  expected,  the  NAVSPACECOM  code  returns  the  SING  =  0 
error  and  does  not  continue  further.  The  SVD  algorithm  on  the  other  hand  returns. 
Further,  upon  analysis  of  the  factored  structure,  the  reconstructed  matrices  are  within  one 
order  of  magnitude  of  the  original  input.  The  differences  are  due  to  roundoff  error. 

While  the  results  presented  in  the  output  of  the  SVD  algorithm  may  not  look  so 
appealing  at  first  glance,  it  is  important  to  note, that  the  present  Gauss-Jordan, 
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NAVSPACECOM  routines  return  no  output  at  all.  This  is  much  more  encouraging.  As  a 
minimum,  even  with  the  ridiculously  singular  nearly  impossible  input  matrices,  the  SVD 
routine  successfully  returned  the  best  approximation  as  the  desired  solution.  At  this  point 
the  differential  correction  subroutine  now  possesses  an  excellent  starting  point  to  continue 
processing  this  observational  input  set.  The  nonlinear  method  may  then  still  converge. 

Had  the  NAVSPACECOM  routine  been  allowed  to  process  this  set,  the  entire  set  would 
have  been  discarded  based  solely  on  the  coincidental  fact  that  the  input  matrix  derivation  is 
singular  beyond  the  computer’s  ability  to  distinguish  otherwise.  The  macro  effect  of  this 
NAVSPACECOM  code  shortfall  is  that  potentially  accurately  derived  observational  data 
is  now  discarded.  This  will  certainly  effect  the  long-term  distribution  of  the  accumulated 
error  in  the  orbital  element  set. 

B.  RECOMMENDATIONS 

The  potential  for  improvement  in  NAVSPACECOMs  differential  correction 
process  through  integration  of  an  SVD  algorithm  which  replace  existing  subroutines 
AMA03,  AMA04,  and  AMA06  exists.  The  only  way  to  verify  the  magnitude  of 
improvement,  however,  is  to  implement  a  test  routine.  The  evidence  presented  herein 
indicates  that  further  study  should  be  pursued. 

C.  CONCLUSION 

The  results  of  this  study  are  not  all  inclusive.  Only  the  worst  possible  input  matrix 
conditions  were  modeled.  The  frequency  of  occurrence  of  this  type  of  input  would 


certainly  have  to  be  taken  into  consideration  when  deciding  whether  to  upgrade  the 
present  process.  Differential  Correction  involves  highly  complex  series  of  numerical 
routines  each  with  its  own  influence  on  the  over  all  solution  process.  My  findings  indicate 
that  some  improvement  may  occur  through  reduced  processing  time  and  the  numerical 
error  of  the  actual  values  returned  in  the  update  state  vectors.  Additionally,  it  is  worth 
pointing  out  that  many  modem  statistical  regression-fitting  packages  have  begun  to  use 
more  sophisticated  linear  algebra  solvers  for  similar  reasons  as  those  addressed  herein. 

While  at  best  the  SVD  algorithm  is  approximately  24  times  more  computationally 
expensive  in  flops,  it  has  the  advantage  of  not  requiring  cyclic  preconditioning  to  start  the 
solution  process.  The  present  NAVSPACECOM  routine  has  the  potential  to  cycle  for 
significant  periods  prior  to  achieving  SING  =  0  stop  criteria.  If  significant  numbers  of 
orbital  element  sets  cause  the  AMA06  routine  to  process  several  times  through  before 
either  failing  or  calling  the  solver,  the  relative  difference  in  run  time  between  the  two 
methods  may  be  very  small.  This  aspect  of  the  study  will  require  actually  run  time 
information  to  determine  the  flop  count  delta. 

SVD  solutions  exhibit  less  error  in  each  iterate  as  a  result  of  their  “ best 
approximation  It  is  possible  that  over  time,  the  resulting  thorough  evaluation  of  the 
observational  data  by  an  SVD  synthesizing  algorithm  could  effectively  decrease  the  error 
in  actual  known  position  information  as  measured  by  the  difference  in  the  orbital  element 
set  and  the  corresponding  test  data. 
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APPENDIX  A 


Test  Case  1.  Singular  A  Matrix 

Highly  singular  A  matrix:  col(4)  =  (col(l)+col(3))/2  and  col(8)  =  (col(2)+col(7))/2 

Constructed  random  uniform  (0,1) 


0.9688 

0.1310 

0.5620 

Matrix  A: 

0.7654  0.5979 

0.0631 

0.7666 

0.4488 

0.3557 

0.9408 

0.3193 

0.3375 

0.9492 

0.2642 

0.6661 

0.8035 

0.0490 

0.7019 

0.3749 

0.2120 

0.2888 

0.9995 

0.1309 

0.4164 

0.7553 

0.8477 

0.8678 

0.8116 

0.8888 

0.2120 

0.0954 

0.4715 

0.8948 

0.2093 

0.3722 

0.6335 

0.1016 

0.4984 

0.0149 

0.1121 

0.2861 

0.4551 

0.0737 

0.1799 

0.0653 

0.2905 

0.2882 

0.3716 

0.2512 

0.0811 

0.1998 

0.2255 

0.2343 

0.6728 

0.8167 

0.4489 

0.9327 

0.8511 

0.0495 

0.4911 

0.9331 

0.9580 

0.9855 

0.9183 

3.4538 

Symmetric,  nxn  matrix  P 

2.2680  1.7824  2.6181  2.6412 

=ata 

1.9559 

• 

2.2782 

2.2731 

2.2680 

3.0953 

1.5425 

1.9052 

2.7916 

2.2445 

1.9391 

2.5172 

1.7824 

1.5425 

1.4978 

1.6401 

1.6543 

1.0673 

1.0142 

1.2783 
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2.6181 

1.9052 

1.6401 

2.1291 

2.1478 

1.5116 

1.6462 

1.7757 

2.6412 

2.7916 

1.6543 

2.1478 

3.0720 

1.8867 

2.3445 

2.5680 

1.9559 

2.2445 

1.0673 

1.5116 

1.8867 

2.8209 

1.9602 

2.1023 

2:2782 

1.9391 

1.0142 

1.6462 

2.3445 

1.9602 

2.7791 

2.3591 

2.2731 

2.5172 

1.2783 

1.7757 

2.5680 

2.1023 

2.3591 

2.4382 

Decomposition  Matrices: 

Matrix  U 


-0.405972 

0.466532 

-0.366483 

0.180020 

-0.532409 

-0.026839 

-0.006785 

0.408157 

-0.387027 

-0.208141 

0.614596 

-0.035119 

-0.301571 

-0.411911 

0.408982 

0.006625 

-0.238695 

0.389374 

0.258549 

0.159762 

0.718607 

-0.126578 

-0.006750 

0.408139 

-0.322333 

0.427960 

-0.053982 

0.169871 

0.093150 

-0.075934 

0.013911 

-0.816422 

-0.404391 

0.039450 

0.240801 

-0.441264 

-0.010121 

0.762960 

-0.001406 

0.000463 

-0.326910 

-0.511282 

-0.089709 

0.724015 

0.092756 

0.301448 

-0.000586 

0.000171 

-0.346071 

-0.284110 

-0.593342 

-0.385717 

0.303011 

-0.199784 

0.408545 

0.006797 

-0.366547 

-0.246130 

0.010611 

-0.210450 

0.000580 

-0.308105 

-0.815805 

-0.013969 

Diagonal  of  Matrix  2 

16.942978 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

Matrix  V-Transpose 
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-0.405972 

-0.387027 

-0.238695 

-0.322333 

-0.404391 

-0.326910 

-0.346070 

-0.366547 

0.467637 

-0.210013 

0.388574 

0.428113 

0.038740 

-0.511046 

-0.282277 

-0.246149 

-0.370840 

0.614798 

0.254242 

-0.058314 

0.255496 

-0.115354 

-0.581027 

0.016870 

0.169917 

-0.013430 

0.165458 

0.167672 

-0.432943 

0.720162 

-0.406542 

-0.210019 

-0.575737 

-0.269208 

0.669202 

0.201179 

-0.065735 

0.073310 

0.312403 

0.021403 

0.261994 

-0.244425 

0.394398 

-0.727229 

0.369304 

0.160414 

-0.054058 

-0.148875 

-0.223596 

-0.360038 

-0.310865 

0.334620 

0.66448 6 

0.262515 

-0.175497 

-0.265547 

0.000238 

-0.407559 

0.000561 

-0.000351 

-0.001333 

-0.000496 

-0.407884 

0.817022 

Check  product  against  original  matrix: 

Original  P  Matrix: 

3.453800 

2.268000 

1.782400 

2.618100 

2.641200 

1.955900 

2.278200 

2.273100 

2.268000 

3.095300 

1.542500 

1.905200 

2.791600 

2.244500 

1.939100 

2.517200 

1.782400 

1.542500 

1.497800 

1.640100 

1.654300 

1.067300 

1.014200 

1.278300 

2.618100 

1.905200 

1.640100 

2.129100 

2.147800 

1.511600 

1.646200 

1.775700 

2.641200 

2.791600 

1.654300 

2.147800 

3.072000 

1.886700 

2.344500 

2.568000 

1.955900 

2.244500 

1.067300 

1.511600 

1.886700 

2.820900 

1.960200 

2.102300 

2.278200 

1.939100 

1.014200 

1.646200 

2.344500 

1.960200 

2.779100 

2.359100 

2.273100 

2.517200 

1.278300 

1.775700 

2.568000 

2.102300 

2.359100 

2.438200 

Product  Matrix  =  U*£*(V-Transpose): 
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2.792422  2.662113 
2.662113  2.537884 

1.641832  1.565216 

2.217128  2.113665 
2.781547  2.651745 
2.248604  2.143672 
2.380401  2.269319 
2.521245  2.403590 


1.641832  2.217128 
1.565216  2.113665 
0.965332  1.303582 
1.303582  1.760356 

1.635438  2.208493 
1.322089  1.785347 

1.399580  1.889991 

1.482391  2.001818 


2.781547 

2.248604 

2.651745 

2.143672 

1.635438 

1.322089 

2.208493 

1.785347 

2.770714 

2.239847 

2.239847 

1.810694 

2.371130 

1.916823 

2.511426 

2.030238 

2.380400 

2.521245 

2.269318 

2.403590 

1.399580 

1.482391 

1.889991 

2.001818 

2.371130 

2.511426 

1.916823 

2.030238 

2.029173 

2.149235 

2.149235 

2.276402 
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APPENDIX  B 


Test  Case  2.  Extreme  Singularity 

Singular  P  matrix:  col(l)  =  col(2)+col(3)  +  . . .  +  col(8) 


1.000000 

0.900000 

0.090000 

0.009000 

0.000900 

0.000090 

0.000009 

0.000001 

0.900000 

0.900000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.090000 

0.000000 

0.090000 

o.oooooo 

0.000000 

0.000000 

0.000000 

0.000000 

0.009000 

0.000000 

0.000000 

0.009000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000900 

0.000000 

0.000000 

0.000000 

0.000900 

0.000000 

0.000000 

0.000000 

0.000090 

0.000000 

0.000000 

0.000000 

0.000000 

0.000090 

0.000000 

0.000000 

0.000009 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000009 

0.000000 

0.000001 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000001 

Decomposition  Matrices: 

Matrix  U 

-0.726830  -0.298545  0.377237 

0.221507 

-0.175290  0.150237 

-0.127672 

0.348738 

-0.685806 

0.302705  -0.443863 

-0.223180 

0.175290  -0.150237 

0.127672 

-0.348738 

-0.037087 

0.334128  0.811978 

-0.190945 

0.175288  -0.150237 

0.127672 

-0.348738 

-0.003546 

-0.836417  0.031208 

-0.329137 

0.174926  -0.150235 

0.127672 

-0.348738 

-0.000353 

-0.089017  -0.019792 

0.864991 

0.288233  -0.149928 

0.127670 

-0.348738 

45 


-0.000035  -0.008940  -0.002062  0.089682  -0.886388  -0.261431  0.127375  -0.348738 


-0.000004  -0.000894  -0.000207  0.008994  -0.091411  0.899086  0.248332  -0.348622 

0.000000  -0.000099  -0.000023  0.001000  -0.010180  0.102642  -0.916846  -0.385685 

Diagonal  of  Matrix  X 

1.860190  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000 

Matrix  V-Transpose 

-0.676744  -0.731343  -0.084492  -0.003917  -0.000497  -0.000049  -0.000005  -0.000001 

0.634569  -0.521247  -0.570616  -0.004328  -0.001723  -0.000168  -0.000017  -0.000002 

0.187067  -0.224736  0.406687  0.865145  0.024989  0.002324  0.000231  0.000026 

0.322910  -0.377922  0.708234  -0.499791  -0.038919  -0.001936  -0.000172  -0.000019 

0.008663  -0.010371  0.016396  -0.041196  0.998458  0.030457  0.002268  0.000243 

0.000000  0.000017  0.000174  0.001730  0.030609  -0.999175  -0.026540  -0.002122 

0.000000  -0.000001  -0.000015  -0.000147  -0.001470  -0.026647  0.999247  0.028168 

0.000000  0.000000  0.000001  0.000014  0.000138  0.001378  0.028215  -0.999601 


Check  product  against  original  matrix: 
Original  Matrix  P: 


1.000000 

0.900000 

0.090000 

0.009000 

0.000900 

0.000090 

0.000009 

0.000001 

0.900000 

0.900000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.090000 

0.000000 

0.090000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.009000 

0.000000 

0.000000 

0.009000 

0.000000 

0.000000 

0.000000 

0.000000 
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0.000900 

0.000000 

0.000000 

0.000000 

0.000900 

0.000000 

0.000000 

0.000000 

0.000090 

0.000000 

0.000000 

0.000000 

0.000000 

0.000090 

0.000000 

0.000000 

0.000009 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000009 

0.000000 

0.000001 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000001 

Product  U: 

*Z*(V-Transpose): 

0.914987 

0.988808 

0.114237 

0.005295 

0.000672 

0.000067 

0.000007 

0.000001 

0.863342 

0.932996 

0.107789 

0.004997 

0.000634 

0.000063 

0.000006 

0.000001 

0.046687 

0.050454 

0.005829 

0.000270 

0.000034 

0.000003 

0.000000 

0.000000 

0.004464 

0.004824 

0.000557 

0.000026 

0.000003 

0.000000 

0.000000 

0.000000 

0.000444 

0.000480 

0.000055 

0.000003 

0.000000 

0.000000 

0.000000 

0.000000 

0.000044 

0.000048 

0.000006 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000004 

0.000005 

0.000001 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000001 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 
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